
#ifndef ROPUFU_MULTI_SENSOR_MONTE_CARLO_ENGINE_HPP_INCLUDED
#define ROPUFU_MULTI_SENSOR_MONTE_CARLO_ENGINE_HPP_INCLUDED

#include "../information/affected_range.hpp"
#include "../information/information_none.hpp"
#include "../information/information_exact.hpp"
#include "../information/information_at_most.hpp"
#include "../more_math.hpp"
#include "dumper.hpp"
#include "loop.hpp"

#include <aftermath/probability/dist_normal.hpp>
#include <aftermath/probability/empirical_measure.hpp>
#include <aftermath/random/sampler_normal.hpp>
#include <aftermath/template_math.hpp>

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <map>
#include <random>
#include <string>
#include <vector>

namespace ropufu
{
    namespace monte_carlo
    {
        namespace detail
        {
            template <typename t_distribution_type_type>
            struct sampler_switch
            {
                typedef t_distribution_type_type type;
            };

            template <>
            struct sampler_switch<aftermath::probability::dist_normal>
            {
                typedef aftermath::random::sampler_normal_from_engine<1024, std::default_random_engine>::type type;
            };
        }

        template <typename t_statistic_type, typename t_distribution_type>
        struct engine
        {
        };

        template <typename t_statistic_type>
        struct engine<t_statistic_type, aftermath::probability::dist_normal>
        {
            typedef t_statistic_type statistic_type;
            //typedef typename statistic_type::information_type info_type;
            //typedef typename statistic_type::window_type window_type;

            typedef aftermath::probability::dist_normal distribution_type;
            typedef engine<statistic_type, distribution_type> type;

            static const std::size_t count_streams = statistic_type::count_streams;

            //typedef std::int_fast32_t stopping_time_type;
            //typedef std::int_least32_t stopping_time_type;
            typedef double stopping_time_type;
            typedef std::default_random_engine uniform_engine_type;
            typedef dumper<stopping_time_type> dumper_type;
            typedef detail::loop<statistic_type, stopping_time_type> loop_type;

            typedef aftermath::probability::dist_normal loglikelihood_dist_type;
            //typedef std::normal_distribution<double> loglikelihood_dist_type;

            typedef detail::sampler_switch<loglikelihood_dist_type>::type sampler_type;
            typedef aftermath::probability::empirical_measure<stopping_time_type, double> empirical_time_type;
            typedef aftermath::probability::empirical_measure<double, double> empirical_double_type;
            typedef aftermath::probability::empirical_measure<bool> empirical_bool_type;

        private:
            double m_mu, m_theta, m_eta;
            uniform_engine_type m_uniform_engine;
            dumper_type m_dumper;
            sampler_type m_sampler_before;
            sampler_type m_sampler_after;

            //static std::uint_fast32_t make_seed(std::uint_fast32_t seed_override)
            //{
            //    std::uint_fast32_t seed = static_cast<std::uint_fast32_t>(std::chrono::high_resolution_clock::now().time_since_epoch().count());
            //    if (seed_override != 0) seed = seed_override;
            //    return seed;
            //}

            void make_seed(std::uint_fast32_t seed_override)
            {
                if (seed_override != 0) this->m_uniform_engine = uniform_engine_type(seed_override);
                else
                {
                    //std::array<typename std::random_device::result_type, uniform_engine_type::state_size> seed_data;
                	std::array<typename std::random_device::result_type, 4> seed_data;
                    std::random_device r;
                    std::generate_n(seed_data.data(), seed_data.size(), std::ref(r));
                    std::seed_seq seq(std::begin(seed_data), std::end(seed_data));

                    this->m_uniform_engine = uniform_engine_type(seq);
                }
            }

        public:
            explicit engine(double mu, double theta, double eta, std::uint_fast32_t seed_override = 0)
                : m_mu(mu), m_theta(theta), m_eta(eta), m_dumper(mu, theta, eta)
            {
                this->make_seed(seed_override);

                // Unaffected sensors.
                double llr_mu_before = (mu - theta) * (theta - mu) / 2;
                double llr_sigma_before = (theta - mu > 0.0) ? (theta - mu) : (mu - theta);
                auto   llr_dist_before = loglikelihood_dist_type(llr_mu_before, llr_sigma_before);
                // Affected sensors.
                double llr_mu_after = (2 * eta - theta - mu) * (theta - mu) / 2;
                double llr_sigma_after = (theta - mu > 0.0) ? (theta - mu) : (mu - theta);
                auto   llr_dist_after = loglikelihood_dist_type(llr_mu_after, llr_sigma_after);

                this->m_sampler_before = sampler_type(llr_dist_before);
                this->m_sampler_after = sampler_type(llr_dist_after);
            }

            // WARNING: <threshold> has to be sorted.
            void execute(std::uint_fast64_t n, statistic_type& procedure,
                std::size_t count_affected, const std::vector<double>& threshold,
                std::vector<double>& expected_time, std::vector<double>& expected_aux,
                bool do_stop_on_reset = false, bool is_quiet = false, bool do_dump = false)
            {
                double expected_shift = this->m_theta - this->m_mu;
                double actual_shift = this->m_eta - this->m_mu;

                std::cout << "~~ K = " << count_streams << ", N = " << n << ", theta = " << expected_shift << ": " << procedure << " ~~" << std::endl;

                std::string filename = this->m_dumper.build_filename(procedure, n, count_affected, threshold, do_stop_on_reset);
                if (do_dump && dumper_type::does_exist(filename)) return; // TODO: load data from file.

                //std::vector<loop_type> threaded_loop(static_cast<std::size_t>(0));
                //threaded_loop.reserve(count_threads);
                //for (std::size_t i = 0; i < count_threads; i++) threaded_loop.emplace_back(n, count_affected, threshold);

                loop_type loop(n, count_affected, threshold);
                loop.execute(procedure, do_stop_on_reset, this->m_uniform_engine, this->m_sampler_before, this->m_sampler_after);

                //std::vector<std::thread> threads(static_cast<std::size_t>(0));
                //threads.reserve(count_threads);
                //
                //for (std::size_t i = 0; i < count_threads; i++)
                //{
                //    loop_type& loop = threaded_loop[i];
                //    // <procedure> is captured by value intentionally.
                //    threads.emplace_back(
                //        [&loop, procedure, do_stop_on_reset, &this->m_uniform_engine, &this->m_sampler_before, &this->m_sampler_after]
                //    {
                //        loop.execute(procedure, do_stop_on_reset, this->m_uniform_engine, this->m_sampler_before, this->m_sampler_after);
                //    });
                //}
                //for (auto& t : threads) t.join();
                //threads.clear();


                const std::map<double, double> z_confidence{
                    { 0.90, 1.644853626951 },
                    { 0.95, 1.959963984540 },
                    { 0.98, 2.326347874041 },
                    { 0.99, 2.575829303549 },
                    { 0.995, 2.807033768344 }
                };

                std::vector<empirical_time_type> run_length = loop.run_length();
                std::vector<empirical_double_type> aux = loop.aux();
                std::vector<empirical_double_type> reset_flag = loop.reset_flag();

                //for (std::size_t i = 1; i < count_threads; i++)
                //{
                //    loop_type& loop = threaded_loop[i];
                //    for (std::size_t j = 0; j < threshold.size(); j++)
                //    {
                //        run_length[j].merge(loop.run_length()[j]);
                //        aux[j].merge(loop.aux()[j]);
                //        reset_flag[j].merge(loop.reset_flag()[j]);
                //    }
                //}

                std::vector<double> err_time(threshold.size());
                for (std::size_t j = 0; j < threshold.size(); j++)
                {
                    expected_time[j] = run_length[j].get_mean();
                    expected_aux[j] = aux[j].get_mean();

                    err_time[j] = z_confidence.at(0.95) * run_length[j].compute_standard_deviation() / std::sqrt(static_cast<double>(n));
                }

                if (!is_quiet)
                {
                    for (std::size_t j = 0; j < threshold.size(); j++)
                    {
                        std::cout << "a = " << threshold[j] << ": ";

                        if (actual_shift == 0.0 || count_affected == 0) std::cout << "ARL";
                        else if (actual_shift == expected_shift) std::cout << "ADD ^ " << count_affected;
                        else std::cout << "True shift of " << actual_shift << " delay ^ " << count_affected;

                        std::cout << " = " << expected_time[j] << " \\pm " << err_time[j] << std::endl;
                        if (do_stop_on_reset)
                        {
                            std::cout << "AUX = " << expected_aux[j] << std::endl;
                            std::cout << "P(stop on reset) = " << reset_flag[j].get_mean() << std::endl;
                        }
                    }
                }

                if (do_dump) this->m_dumper.dump(filename, n, count_affected, threshold, run_length, aux, do_stop_on_reset);
            }

            template <typename t_information_type>
            void change_of_measure(statistic_type& proc, double theta, std::vector<double>& threshold, std::uint_fast64_t n,
                bool do_stop_on_reset = true, bool is_quiet = false, bool do_dump = false)
            {
                typedef t_information_type info_type;
                std::vector<double> temp(threshold.size());
                std::vector<double> aux(threshold.size());

                std::vector<double> expected_time_infinity(threshold.size());
                std::vector<double> expected_aux_zero(threshold.size());
                for (std::size_t j = 0; j < threshold.size(); j++) expected_aux_zero[j] = 0.0;

                // ~~ ARL ~~
                this->execute(n, proc, 0, threshold, expected_time_infinity, temp, do_stop_on_reset, is_quiet, do_dump); // ARL
                // ~~ P(cross up) ~~
                double denom = 0.0;
                for (std::size_t k = information::affected_range<info_type>::value_from; k <= information::affected_range<info_type>::value_to; k++)
                {
                    std::size_t weight = more_math::n_choose_k(count_streams, k);
                    if (!is_quiet) std::cout << "---- (" << k << ") weight = " << weight << std::endl;

                    this->execute(n, proc, k, threshold, temp, aux, do_stop_on_reset, false, do_dump); // ADD _0 ^k
                    for (std::size_t j = 0; j < threshold.size(); j++) expected_aux_zero[j] += weight * aux[j];
                    denom += (weight * k) / count_streams;
                }
                for (std::size_t j = 0; j < threshold.size(); j++) expected_aux_zero[j] /= denom;
                if (!is_quiet) std::cout << "---- weight norm = " << denom << std::endl;

                if (!is_quiet)
                {
                    for (std::size_t j = 0; j < threshold.size(); j++)
                    {
                        std::cout << "E _inf (T) = " << expected_time_infinity[j] << std::endl;
                        std::cout << "P _inf (cross up) = " << expected_aux_zero[j] << std::endl;
                        std::cout << "E _inf (S) = " << expected_time_infinity[j] / expected_aux_zero[j] << std::endl;
                    }
                }
            }
        };
    }
}

#endif // ROPUFU_MULTI_SENSOR_MONTE_CARLO_ENGINE_HPP_INCLUDED
